function [gbestX,gbestfitness,gbesthistory]=CS(popsize,dimension,xmax,xmin,~,~,maxiter,Func,FuncId)


gbestX=rand(1,dimension);
Pa = 0.25 ; % Probability of building a new nest(After host bird find exotic bird eggs)
FEs=0;
MaxFEs=10000*dimension;
gbesthistory=[];
gbestfitness=inf;
fitness=inf(1,popsize);
x = rand(popsize,dimension)*(xmax-xmin)+xmin ;  % 初始化寄主的鸟巢Random initial solutions
while 1
    levy_nestPop =  func_levy(x,xmax,xmin) ; % 通过levy飞行产生一个解Generate new solutions by Levy flights


    [gbestX,FEs,gbestfitness,gbesthistory,x,fitness] = func_bestNestPop(gbestX,FEs,gbestfitness,gbesthistory,x,levy_nestPop,Func,FuncId,fitness);  % 与上一代比较，更新适应度较优的鸟巢Choose a best nest among  new and old nests
    rand_nestPop = func_newBuildNest(x,Pa,xmax,xmin); % 根据发现概率舍弃一个鸟巢并建立一个新鸟巢Abandon(Pa) worse nests and build new nests by (Preference random walk )
   
    
    
    [gbestX,FEs,gbestfitness,gbesthistory,x,fitness] = func_bestNestPop(gbestX,FEs,gbestfitness,gbesthistory,x,rand_nestPop,Func,FuncId,fitness) ; %列出当前最佳的鸟巢 Choose a best nest among  new and old nests
    if FEs>=MaxFEs
        break;
    end
end

if FEs<MaxFEs
    gbesthistory(FEs+1:MaxFEs)=gbestfitness;
else
    if FEs>MaxFEs
        gbesthistory(MaxFEs+1:end)=[];
    end
end
end


%% Levy飞行 func_levy.m
 function [ result ] = func_levy( x,xmax,xmin)
%FUNC_LEVY : Update position of nest by using Levy flights
%@author : zhaoyuqiang 
[popsize,dimension] = size(x) ;
% Levy flights by Mantegna's algorithm
beta = 1.5 ;
alpha = 1 ;
sigma_u = (gamma(1+beta)*sin(pi*beta/2)/(beta*gamma((1+beta)/2)*2^((beta-1)/2)))^(1/beta) ;
sigma_v = 1 ;
u = normrnd(0,sigma_u,popsize,dimension) ;%(第一个参数代表均值，sigma参数代表标准差),生成N×D形式的正态分布的随机数矩阵。
v = normrnd(0,sigma_v,popsize,dimension) ;
step = u./(abs(v).^(1/beta)) ;
% alpha = 0.1.*(x(randperm(popsize),:)-x(randperm(popsize),:)); % Bad effect
x = x+alpha.*step ;
% Deal with bounds
x(find(x>xmax)) = xmax ; %#ok<*FNDSB>%查找大于Xmax的元素
x(find(x<xmin)) = xmin ;
result = x ; 
 end


%% 与上一代比较，返回较优的鸟巢 func_bestNestPop.m
function [gbestX,FEs,gbestfitness,gbesthistory,x,fitness] = func_bestNestPop(gbestX,FEs,gbestfitness,gbesthistory,x,new_nestPop,Func,FuncId,fitness)

[gbestX,FEs,gbestfitness,gbesthistory,newfitness]=func_fitness(gbestX,FEs,gbestfitness,gbesthistory,new_nestPop,Func,FuncId);
index = find(fitness>newfitness); %与上一代比较适应度，选择一个适应度更大的更新鸟窝
x(index,:) = new_nestPop(index,:) ;
fitness(index)=newfitness(index);
end

%% 根据发现概率，舍弃一个鸟巢并建立一个新鸟巢 func_newBuildNest.m
function [ x ] = func_newBuildNest( x ,Pa ,xmax,xmin)
%FUNC_NEWBUILDNEST new solutions are generated by using the similarity 
% between the existing eggs/solutions and the host eggs/solutions with a discovery rate pa .
%@author zhaoyuqiang
[popsize,dimension] = size(x) ;
%%根据发现概率发现鸟蛋，舍弃鸟窝
x = x+rand.*heaviside(rand(popsize,dimension)-Pa).*(x(randperm(popsize),:)-x(randperm(popsize),:));
% Deal with bounds
x(find(x>xmax)) = xmax ; %#ok<*FNDSB>建立新的鸟窝
x(find(x<xmin)) = xmin ;
end

%% 计算适应度函数
function [gbestX,FEs,gbestfitness,gbesthistory,objValue] = func_fitness(gbestX,FEs,gbestfitness,gbesthistory,pop,Func,FuncId)

Fitness = Func;
[row,col]=size(pop);
for i=1:row
objValue(i) =  Fitness(pop(i,:)',FuncId);
FEs=FEs+1;
if gbestfitness>objValue(i)
    gbestfitness=objValue(i);
    gbestX=pop(i,:);
end
    gbesthistory(FEs)=gbestfitness;
    fprintf("CS算法,第%d代，最佳适应度 = %e\n",FEs,gbestfitness);
end

end

